Loading libraries

library(dada2); packageVersion("dada2")
## [1] '1.18.0'
library(ggplot2); packageVersion("ggplot2")
## [1] '3.3.3'
library(phyloseq); packageVersion("phyloseq")
## [1] '1.34.0'
library(data.table); packageVersion("data.table")
## [1] '1.13.6'
library(dplyr)
library(plyr)
library(magrittr)
library(RColorBrewer)
library(picante)
library(igraph)
library(cowplot)
library(plotly)
library(lemon)
library(GGally)
library(randomcoloR)
library(knitr)
library(kableExtra)
library(ggpmisc)
library(ggpubr)
library(decontam)
library(stringr)
library(tibble)

Importing data

SVt = readRDS("../seqtab_final.rds")
# SVt_cnn = collapseNoMismatch(SVt, minOverlap = 20,
#                               identicalOnly = FALSE, vec = TRUE, verbose = TRUE)
# saveRDS(SVt_cnn, "../seqtab_final_cnn.rds")
# SVt = readRDS("../seqtab_final_cnn.rds")
SVt = SVt[!rownames(SVt) %in% "Undetermined", ]

Tax = readRDS("../tax_final.rds")

metadata<-read.csv("../../220319_for_graph.csv", header=TRUE)
# metadata[, 10:47]  <- as.numeric(metadata[, 10:47])
metadata$Snowpack_Layer = gsub("Glacial Ice", "GL ICE", metadata$Snowpack_Layer)
metadata$Snowpack_Layer = gsub("Sup. Ice", "SUP ICE", metadata$Snowpack_Layer)
metadata$Snowpack_Layer = gsub("Top", "TOP", metadata$Snowpack_Layer)
metadata$Snowpack_Layer = gsub("Mid", "MID", metadata$Snowpack_Layer)
metadata$Snowpack_Layer = gsub("Basal", "BASAL", metadata$Snowpack_Layer)

order<-rownames(SVt)
metadata$Sample_ID=as.character(metadata$Sample_ID)
metadata_ord<-metadata[match(order, metadata$Sample_ID),]
rownames(metadata_ord) <- metadata_ord$Sample_ID
metadata_ord$Sample_Name <- factor(metadata$Sample_Name, levels=unique(metadata$Sample_Name))
ps <- phyloseq(otu_table(SVt, taxa_are_rows=FALSE),
               sample_data(metadata_ord),
               tax_table(Tax))
ps.b <- prune_taxa(taxa_sums(ps) > 0, ps)
ps.b = subset_taxa(ps.b, Kingdom %in% c("Archaea", "Bacteria"))
ps.b = subset_samples(ps.b, Fox_Aspect != "N_NE_D")

Testing batch effects

#remove samples with 0 reads
ps.batch = prune_samples(sample_sums(ps.b)>=1, ps.b)

readsumsdf_samples<-data.frame(nreads = sample_sums(ps.batch),
                               samples = rownames(as.data.frame(sample_sums(ps.batch))),
                               batch = sample_data(ps.batch)$Batch.no.)
readsumsdf_samples_f = readsumsdf_samples[- grep("Control", readsumsdf_samples$batch),]

my_comparisons=list(c("1","2"),c("1","3"),c("1","4"),c("1","5"),
                    c("2","3"), c("2","4"), c("2","5"),
                    c("3","4"), c("3","5"),
                    c("4","5"))

ggplot(readsumsdf_samples_f, aes(x=batch, y=nreads,
                                 color=batch, fill=batch)) +
  geom_violin(alpha=0.3) +
  geom_jitter(shape=16, position = "jitter") +
  stat_compare_means(comparisons = my_comparisons,
                     method = "t.test",
                     symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1), 
                       symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
  labs(color="Batch Number") +
  guides(fill=FALSE) +
  xlab("Extraction batch") +
  ylab("Number of reads obtained") +
  theme(axis.text.x = element_text(vjust = 1,
                                   hjust = 1,
                                   size = 12),
        strip.text = element_text(size = 12, margin = margin()))

1 - Library sizes

sort(sample_sums(ps.b))
##     C05-cDNA-29     D12-cDNA-MM     C10-cDNA-34 A02-cDNA-MM-H20     C08-cDNA-32 
##               0               0               5               7              25 
##     C07-cDNA-31     B07-cDNA-17 E01-cDNA-MM-H20     B11-cDNA-22     B09-cDNA-19 
##              27              43              55             123             187 
##     D01-cDNA-37     A01-cDNA-MM      A11-cDNA-9     D10-cDNA-46     B10-cDNA-21 
##             262             290            2320            3844            9170 
##      A10-cDNA-8     B05-cDNA-15       B08-DNA-4       B11-DNA-7     B08-cDNA-18 
##           15595           15781           16781           29350           29808 
##     C03-cDNA-27      A08-cDNA-6      A06-cDNA-4      A09-cDNA-7      A07-cDNA-5 
##           31806           33041           37485           40604           44392 
##       A11-DNA-9      A05-cDNA-3       A08-DNA-6     B02-cDNA-12       D08-DNA-8 
##           50662           52562           53872           57398           57493 
##      E02-DNA-MM     C02-cDNA-26     D07-cDNA-43       D06-DNA-6       D02-DNA-2 
##           58537           59570           59805           60653           61230 
##      A01-DNA-MM       B12-DNA-8     B12-cDNA-23      A12-DNA-10       C07-DNA-5 
##           61807           62146           62235           63891           64148 
##     B04-cDNA-14       A07-DNA-5  A02-DNA-MM-H20      A04-cDNA-2     B03-cDNA-13 
##           68029           69142           69499           70613           76895 
##       D09-DNA-9      B02-DNA-12       A05-DNA-3     C04-cDNA-28       C05-DNA-3 
##           80028           83134           84732           86983           87568 
##       C04-DNA-2     D06-cDNA-42       A09-DNA-7     B01-cDNA-11      B04-DNA-14 
##           96798          101335          108915          110464          113222 
##    C03-DNA-B2-1      B03-DNA-13       C08-DNA-6     A12-cDNA-10     B06-cDNA-16 
##          114048          118548          122463          123046          123687 
##     C12-cDNA-36     D02-cDNA-38       D03-DNA-3       C06-DNA-4  E03-DNA-MM-H20 
##          123719          137347          138813          139808          140502 
##       D07-DNA-7      C12-DNA-10       D04-DNA-4       B09-DNA-5     D08-cDNA-44 
##          140668          148596          151009          155211          158902 
##       A06-DNA-4       C01-DNA-9       B10-DNA-6       A10-DNA-8       C09-DNA-7 
##          166098          169126          170463          170498          173577 
##      B01-DNA-11       E01-DNA-4     C06-cDNA-30     C09-cDNA-33     D03-cDNA-39 
##          174897          175874          188016          190708          192008 
##      A03-cDNA-1       A03-DNA-1       B06-DNA-2     C11-cDNA-35       A04-DNA-2 
##          193167          196436          198649          198996          200229 
##       B07-DNA-3     D04-cDNA-40       C10-DNA-8     D11-cDNA-47    D01-DNA-B4-1 
##          204577          204668          212618          220186          238199 
##     D09-cDNA-45       D11-DNA-2       D05-DNA-5    D10-DNA-B5-1     C01-cDNA-25 
##          240324          252560          267031          280585          293592 
##     D05-cDNA-41       D12-DNA-3 
##          302761          312011
df <- as.data.frame(sample_data(ps.b)) # Put sample_data into a ggplot-friendly data.frame
df$LibrarySize <- sample_sums(ps.b)
df <- df[order(df$LibrarySize),]
df$Index <- seq(nrow(df))
ggplot(data=df, aes(x=Index, y=LibrarySize, color=Sample_type)) + geom_point()

2 - Evaluating dataset contamination with decontam

#remove samples with 0 reads
ps.b = prune_samples(sample_sums(ps.b)>=1, ps.b)

sample_data(ps.b)$is.neg <- sample_data(ps.b)$Sample_type == "CONTROL"
contamdf0.5.prev <- isContaminant(ps.b, method="prevalence", threshold = 0.1,
                               neg="is.neg")

ps.final.noncontam <- prune_taxa(!contamdf0.5.prev$contaminant, ps.b)
ps.final.noncontam.noconts = subset_samples(ps.final.noncontam, Sample_type != "CONTROL")

ps.final.noncontam.noconts.nodoubs = subset_samples(ps.final.noncontam.noconts,
    Sample_ID != "D02-DNA-2" & Sample_ID != "B11-DNA-7" & Sample_ID != "C10-cDNA-34" & Sample_ID != "A08-cDNA-6")

ps.final.noncontam.noconts.nodoubs = prune_samples(sample_sums(ps.final.noncontam.noconts.nodoubs)>=1000,
                                                   ps.final.noncontam.noconts.nodoubs)

ps.final <- prune_taxa(taxa_sums(ps.final.noncontam.noconts.nodoubs) > 0, ps.final.noncontam.noconts.nodoubs)
ps.final = prune_samples(sample_sums(ps.final)>=1, ps.final)

ps.final.r <-  transform_sample_counts(ps.final, function(x) x / sum(x))
ps.final.r.f = filter_taxa(ps.final.r, function(x) sum(x) > .001, TRUE)
ps.final.r
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 35534 taxa and 65 samples ]
## sample_data() Sample Data:       [ 65 samples by 47 sample variables ]
## tax_table()   Taxonomy Table:    [ 35534 taxa by 6 taxonomic ranks ]
ps.final.r.f
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 2886 taxa and 65 samples ]
## sample_data() Sample Data:       [ 65 samples by 47 sample variables ]
## tax_table()   Taxonomy Table:    [ 2886 taxa by 6 taxonomic ranks ]
ps.final
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 35534 taxa and 65 samples ]
## sample_data() Sample Data:       [ 65 samples by 47 sample variables ]
## tax_table()   Taxonomy Table:    [ 35534 taxa by 6 taxonomic ranks ]
sum(rowSums(otu_table(ps.final)))
## [1] 6750655
gl_ice_spls = as.character(unique(sample_data(ps.final)[str_detect(sample_data(ps.final)$Sample_Name, "Gl ice"),]$Sample_Name))

Setting order for all plots

sample_data(ps.final.r)$Month = factor(sample_data(ps.final.r)$Month, 
                                 levels = c("April","June","Early July", "Late July"))
sample_data(ps.final.r)$Fox_Aspect = factor(sample_data(ps.final.r)$Fox_Aspect, 
                                 levels = c("NW_AWS","N_NE","SWS1SE"))
sample_data(ps.final.r)$Snowpack_Layer = factor(sample_data(ps.final.r)$Snowpack_Layer, 
                                 levels = c("TOP","MID","BASAL","SUP ICE","GL ICE"))

sample_data(ps.final)$Month = factor(sample_data(ps.final)$Month, 
                                 levels = c("April","June","Early July", "Late July"))
sample_data(ps.final)$Fox_Aspect = factor(sample_data(ps.final)$Fox_Aspect, 
                                 levels = c("NW_AWS","N_NE","SWS1SE"))
sample_data(ps.final)$Snowpack_Layer = factor(sample_data(ps.final)$Snowpack_Layer, 
                                 levels = c("TOP","MID","BASAL","SUP ICE","GL ICE"))

sample_data(ps.final.r.f)$Month = factor(sample_data(ps.final.r.f)$Month, 
                                 levels = c("April","June","Early July", "Late July"))
sample_data(ps.final.r.f)$Fox_Aspect = factor(sample_data(ps.final.r.f)$Fox_Aspect, 
                                 levels = c("NW_AWS","N_NE","SWS1SE"))
sample_data(ps.final.r.f)$Snowpack_Layer = factor(sample_data(ps.final.r.f)$Snowpack_Layer, 
                                 levels = c("TOP","MID","BASAL","SUP ICE","GL ICE"))

3a - Phyla- and Proteobacterial class-level stacked barchart by DNA_or_cDNA~Fox_Aspect, “GL ICE” removed, “Gl ice” samples removed, samples “30 July_SWS1SE_Slush_R”, “30 July_SWS1SE_Slush_D” also removed

ps.final.r.glom <- tax_glom(subset_samples(ps.final.r, 
                             Snowpack_Layer != "GL ICE" & 
                             Sample_Name != "30 July_SWS1SE_Slush_R" &
                             Sample_Name != "30 July_SWS1SE_Slush_D" &
                             !(Sample_Name %in% gl_ice_spls)), taxrank = 'Phylum', NArm = FALSE)
ps.final.r.glom.psdf <- data.table(psmelt(ps.final.r.glom))
ps.final.r.glom.psdf$Phylum <- as.character(ps.final.r.glom.psdf$Phylum)

ps.final.r.glom.psdf[, mean := mean(Abundance, na.rm = TRUE), 
                 by = "Phylum"]
ps.final.r.glom.psdf[(mean <= 0.01), Phylum := "Other"]
ps.final.r.glom.psdf$Phylum[is.na(ps.final.r.glom.psdf$Phylum) & ps.final.r.glom.psdf$Kingdom=="Bacteria"] <- "Unc. Bacteria"
ps.final.r.glom.psdf$Phylum[is.na(ps.final.r.glom.psdf$Phylum) & ps.final.r.glom.psdf$Kingdom=="Archaea"] <- "Unc. Archaea"
#Removing Proteobacteria from previous table
ps.final.r.glom.psdf.noProteo<-ps.final.r.glom.psdf[!grepl("Proteobacteria", ps.final.r.glom.psdf$Phylum),]

#New table, with Proteobacterial classes only
ps.final.r.ProtCl = subset_taxa(ps.final.r, Phylum == "Proteobacteria")
ps.final.r.ProtClglom <- tax_glom(ps.final.r.ProtCl, taxrank = 'Class', NArm = FALSE)
ps.final.r.ProtClglom.psdf <- data.table(psmelt(ps.final.r.ProtClglom))
ps.final.r.ProtClglom.psdf$Class <- as.character(ps.final.r.ProtClglom.psdf$Class)

ps.final.r.ProtClglom.psdf[, mean := mean(Abundance, na.rm = TRUE), 
                       by = "Class"]
ps.final.r.ProtClglom.psdf[(mean <= 0.001), Class := "Other Proteobacteria"]
ps.final.r.ProtClglom.psdf.noP <- subset(ps.final.r.ProtClglom.psdf, select = -Phylum)

names(ps.final.r.glom.psdf.noProteo)[names(ps.final.r.glom.psdf.noProteo) == 'Phylum'] <- 'Taxa'
names(ps.final.r.ProtClglom.psdf.noP)[names(ps.final.r.ProtClglom.psdf.noP) == 'Class'] <- 'Taxa'
# ps.final.r.glom.psdf.noProteo = subset(ps.final.r.glom.psdf.noProteo, select=-c(Class))
ps.final.r.ProteoClAllPhy <- rbind(ps.final.r.glom.psdf.noProteo, ps.final.r.ProtClglom.psdf.noP)

tol18rainbow=c("#771155", "#AA4488", "#CC99BB", "#114477", "#4477AA", "#77AADD", "#117777", "#44AAAA", "#77CCCC", "#777711", "#AAAA44", "#DDDD77", "#774411", "#AA7744", "#DDAA77", "#771122", "#AA4455")

ps.final.r.ProteoClAllPhy$Taxa = factor(ps.final.r.ProteoClAllPhy$Taxa, 
              levels = c("Alphaproteobacteria", "Deltaproteobacteria",
                         "Gammaproteobacteria",
                         "Other Proteobacteria",
                         "Actinobacteria",
                         "Acidobacteria","Bacteroidetes",
                         "Chloroflexi","Cyanobacteria",
                         "Firmicutes","Gemmatimonadetes","Other"))
ggplot(ps.final.r.ProteoClAllPhy, 
       aes(x = Month, y = Abundance, fill = Taxa)) + 
  facet_grid(DNA_or_cDNA~Fox_Aspect,
             scales = "free",
             space = "free") +
  geom_bar(stat = "identity",
           position = "fill") +
  scale_fill_manual(values = tol18rainbow,
                    na.value = "gray40") +
  scale_y_continuous(expand = c(0,0),
                     labels = scales::percent,
                     breaks = c(0.25,0.5,0.75,1)) +
  theme(axis.title.x = element_blank(),
        axis.text.x = element_text(angle = 35, hjust = 1),
        strip.text.x = element_text(face = "bold",
                                    size = 12),
        strip.background = element_rect(size=1),
        legend.text = element_text(size= 12)) + 
  guides(fill = guide_legend(keywidth = 1, keyheight = 1)) +
  ylab("Relative Abundance\n")

3b - Phyla- and Proteobacterial class-level stacked barchart by DNA_or_cDNA~Snowpack_Layer

ggplot(ps.final.r.ProteoClAllPhy, 
                                   aes(x = Month, y = Abundance, fill = Taxa)) + 
  facet_grid(DNA_or_cDNA~Snowpack_Layer,
             scales = "free",
             space = "free") +
  geom_bar(stat = "identity",
           position = "fill") +
  scale_fill_manual(values = tol18rainbow,
                    na.value = "gray40") +
  scale_y_continuous(expand = c(0,0),
                     labels = scales::percent,
                     breaks = c(0.25,0.5,0.75,1)) +
  theme(axis.title.x = element_blank(),
        axis.text.x = element_text(angle = 35, hjust = 1),
        strip.text.x = element_text(face = "bold",
                                    size = 11),
        strip.background = element_rect(size=1),
        legend.text = element_text(size= 12)) + 
  guides(fill = guide_legend(keywidth = 1, keyheight = 1)) +
  ylab("Relative Abundance\n")

3c - Phyla- and Proteobacterial class-level stacked barchart by DNA_or_cDNA~Month

tol17rainbow=c("#771155", "#AA4488", "#CC99BB", "#114477", "#4477AA", "#77AADD", "#117777", "#44AAAA", "#77CCCC", "#777711", "#AAAA44", "#DDDD77", "#774411", "#AA7744", "#DDAA77", "#AA4455")

ggplot(ps.final.r.ProteoClAllPhy,
                                   aes(x = Sample_Name, y = Abundance, fill = Taxa)) +
  facet_grid(DNA_or_cDNA~Month,
             scales = "free",
             space = "free") +
  geom_bar(stat = "identity",
           position = "fill",
           color="black") +
  scale_fill_manual(values = tol17rainbow,
                    na.value = "gray40") +
  scale_y_continuous(expand = c(0,0),
                     labels = scales::percent,
                     breaks = c(0.25,0.5,0.75,1)) +
  theme(axis.title.x = element_blank(),
        axis.text.x = element_text(angle = 35, hjust = 1),
        strip.text.x = element_text(face = "bold",
                                    size = 16),
        strip.background = element_rect(size=1),
        legend.text = element_text(size= 14)) +
  guides(fill = guide_legend(keywidth = 1, keyheight = 1,
                             nrow = )) +
  ylab("Relative Abundance\n")

4 - Top SVs and their sequences

Top20SVs_names = names(sort(taxa_sums(ps.final.r), TRUE)[1:20])
Top20SVs_names_df=as.data.frame(as(tax_table(prune_taxa(Top20SVs_names, ps.final.r)), "matrix")[,2:6])
rownames(Top20SVs_names_df) = c()
Top20SVs_names_df
##            Phylum               Class                 Order            Family
## 1  Proteobacteria Alphaproteobacteria           Rhizobiales  Beijerinckiaceae
## 2  Proteobacteria Gammaproteobacteria Betaproteobacteriales  Burkholderiaceae
## 3  Proteobacteria Alphaproteobacteria           Rhizobiales         Labraceae
## 4  Proteobacteria Alphaproteobacteria      Sphingomonadales Sphingomonadaceae
## 5  Proteobacteria Gammaproteobacteria Betaproteobacteriales  Burkholderiaceae
## 6  Proteobacteria Alphaproteobacteria           Rhizobiales Xanthobacteraceae
## 7  Proteobacteria Alphaproteobacteria      Sphingomonadales Sphingomonadaceae
## 8  Proteobacteria Alphaproteobacteria           Rhizobiales  Beijerinckiaceae
## 9  Proteobacteria Alphaproteobacteria      Sphingomonadales Sphingomonadaceae
## 10 Proteobacteria Alphaproteobacteria       Acetobacterales  Acetobacteraceae
## 11 Actinobacteria      Actinobacteria         Micrococcales Microbacteriaceae
## 12 Proteobacteria Gammaproteobacteria Betaproteobacteriales  Burkholderiaceae
## 13 Proteobacteria Alphaproteobacteria      Sphingomonadales Sphingomonadaceae
## 14 Proteobacteria Gammaproteobacteria Betaproteobacteriales  Burkholderiaceae
## 15 Actinobacteria      Actinobacteria         Micrococcales    Micrococcaceae
## 16 Actinobacteria      Actinobacteria     Corynebacteriales      Nocardiaceae
## 17 Actinobacteria      Actinobacteria            Frankiales              <NA>
## 18 Actinobacteria      Actinobacteria         Micrococcales Cellulomonadaceae
## 19 Actinobacteria      Actinobacteria            Frankiales              <NA>
## 20 Actinobacteria      Actinobacteria            Frankiales              <NA>
##                                         Genus
## 1                                       Bosea
## 2                                 Cupriavidus
## 3                                      Labrys
## 4                                Sphingomonas
## 5  Burkholderia-Caballeronia-Paraburkholderia
## 6                              Bradyrhizobium
## 7                                Sphingomonas
## 8                            Methylobacterium
## 9                                Sphingomonas
## 10                                       <NA>
## 11                            Salinibacterium
## 12                                   Massilia
## 13                               Sphingomonas
## 14                                    Delftia
## 15                               Arthrobacter
## 16                                       <NA>
## 17                                       <NA>
## 18                                       <NA>
## 19                                       <NA>
## 20                                       <NA>
rownames(as.data.frame(as(tax_table(prune_taxa(Top20SVs_names, ps.final.r)), "matrix")[,2:6]))
##  [1] "TGGGGAATATTGGACAATGGGCGCAAGCCTGATCCAGCCATGCCGCGTGAGTGATGAAGGCCTTAGGGTTGTAAAGCTCTTTTGTCCGGGAAGATAATGACTGTACCGGAAGAATAAGCCCCGGCTAACTTCGTGCCAGCAGCCGCGGTAATACGAAGGGGGCTAGCGTTGCTCGGAATCACTGGGCGTAAAGGGCGCGTAGGCGGACTCTTAAGTCGGGGGTGAAAGCCCAGGGCTCAACCCTGGAATTGCCTTCGATACTGGGAGTCTTGAGTTCGGAAGAGGTTGGTGGAACTGCGAGTGTAGAGGTGAAATTCGTAGATATTCGCAAGAACACCGGTGGCGAAGGCGGCCAACTGGTCCGAAACTGACGCTGAGGCGCGAAAGCGTGGGGAGCAAACA"                         
##  [2] "TGGGGAATTTTGGACAATGGGGGCAACCCTGATCCAGCAATGCCGCGTGTGTGAAGAAGGCCTTCGGGTTGTAAAGCACTTTTGTCCGGAAAGAAATCGCGCTGGTTAATACCTGGCGTGGATGACGGTACCGGAAGAATAAGCACCGGCTAACTACGTGCCAGCAGCCGCGGTAATACGTAGGGTGCGAGCGTTAATCGGAATTACTGGGCGTAAAGCGTGCGCAGGCGGTTTTGTAAGACAGGCGTGAAATCCCCGGGCTTAACCTGGGAATTGCGCTTGTGACTGCAAGGCTAGAGTGCGTCAGAGGGGGGTAGAATTCCACGTGTAGCAGTGAAATGCGTAGAGATGTGGAGGAATACCGATGGCGAAGGCAGCCCCCTGGGACGTGACTGACGCTCATGCACGAAAGCGTGGGGAGCAAACA"
##  [3] "TGGGGAATATTGGACAATGGGCGCAAGCCTGATCCAGCCATGCCGCGTGAGTGATGACGGCCTTAGGGTTGTAAAGCTCTTTTAACAGGGACGATAATGACGGTACCTGTAGAATAAGCCCCGGCAAACTTCGTGCCAGCAGCCGCGGTAATACGAAGGGGGCTAGCGTTGTTCGGAATTACTGGGCGTAAAGCGCACGTAGGCGGATTGTTAAGTCGGGGGTGAAATCCTGAGGCTCAACCTCAGAACTGCCTTCGATACTGGCGATCTTGAGTTCGGAAGAGGTTGGTGGAACAGCTAGTGTAGAGGTGAAATTCGTAGATATTAGCTAGAACACCAGTGGCGAAGGCGGCCAACTGGTCCGATACTGACGCTGAGGTGCGAAAGCGTGGGGAGCAAACA"                         
##  [4] "TGGGGAATATTGGACAATGGGCGAAAGCCTGATCCAGCAATGCCGCGTGAGTGATGAAGGCCTTAGGGTTGTAAAGCTCTTTTACCCGGGATGATAATGACAGTACCGGGAGAATAAGCTCCGGCTAACTCCGTGCCAGCAGCCGCGGTAATACGGAGGGAGCTAGCGTTATTCGGAATTACTGGGCGTAAAGCGCACGTAGGCGGCTTTGTAAGTAAGAGGTGAAAGCCTGGTGCTCAACACCAGAACTGCCTTTTAGACTGCATCGCTTGAATCCAGGAGAGGTGAGTGGAATTCCGAGTGTAGAGGTGAAATTCGTAGATATTCGGAAGAACACCAGTGGCGAAGGCGGCTCACTGGACTGGTATTGACGCTGAGGTGCGAAAGCGTGGGGAGCAAACA"                         
##  [5] "TGGGGAATTTTGGACAATGGGGGAAACCCTGATCCAGCAATGCCGCGTGTGTGAAGAAGGCCTTCGGGTTGTAAAGCACTTTTGTCCGGAAAGAAAACCTTTGCACTAATACTGTGAGGGGATGACGGTACCGGAAGAATAAGCACCGGCTAACTACGTGCCAGCAGCCGCGGTAATACGTAGGGTGCGAGCGTTAATCGGAATTACTGGGCGTAAAGCGTGCGCAGGCGGTTCGTTAAGACAGATGTGAAATCCCCGGGCTTAACCTGGGAACTGCATTTGTGACTGGCGAGCTAGAGTATGGCAGAGGGGGGTAGAATTCCACGTGTAGCAGTGAAATGCGTAGAGATGTGGAGGAATACCGATGGCGAAGGCAGCCCCCTGGGCCAATACTGACGCTCATGCACGAAAGCGTGGGGAGCAAACA"
##  [6] "TGGGGAATATTGGACAATGGGCGCAAGCCTGATCCAGCCATGCCGCGTGAGTGATGAAGGCCCTAGGGTTGTAAAGCTCTTTTGTGCGGGAAGATAATGACGGTACCGCAAGAATAAGCCCCGGCTAACTTCGTGCCAGCAGCCGCGGTAATACGAAGGGGGCTAGCGTTGCTCGGAATCACTGGGCGTAAAGGGTGCGTAGGCGGGTCTTTAAGTCAGGGGTGAAATCCTGGAGCTCAACTCCAGAACTGCCTTTGATACTGAAGATCTTGAGTTCGGGAGAGGTGAGTGGAACTGCGAGTGTAGAGGTGAAATTCGTAGATATTCGCAAGAACACCAGTGGCGAAGGCGGCTCACTGGCCCGATACTGACGCTGAGGCACGAAAGCGTGGGGAGCAAACA"                         
##  [7] "TGGGGAATATTGGACAATGGGCGAAAGCCTGATCCAGCAATGCCGCGTGAGTGATGAAGGCCTTAGGGTTGTAAAGCTCTTTTACCCGGGATGATAATGACAGTACCGGGAGAATAAGCTCCGGCTAACTCCGTGCCAGCAGCCGCGGTAATACGGAGGGAGCTAGCGTTATTCGGAATTACTGGGCGTAAAGCGCACGTAGGCGGCTTTGTAAGTAAGAGGTGAAAGCCCAGAGCTCAACTCTGGAATTGCCTTTTAGACTGCATCGCTTGAATCATGGAGAGGTCAGTGGAATTCCGAGTGTAGAGGTGAAATTCGTAGATATTCGGAAGAACACCAGTGGCGAAGGCGGCTGACTGGACATGTATTGACGCTGAGGTGCGAAAGCGTGGGGAGCAAACA"                         
##  [8] "TGGGGAATATTGGACAATGGGCGCAAGCCTGATCCAGCCATGCCGCGTGAGTGATGAAGGCCTTAGGGTTGTAAAGCTCTTTTATCCGGGACGATAATGACGGTACCGGAGGAATAAGCCCCGGCTAACTTCGTGCCAGCAGCCGCGGTAATACGAAGGGGGCTAGCGTTGCTCGGAATCACTGGGCGTAAAGGGCGCGTAGGCGGCGTTTTAAGTCGGGGGTGAAAGCCTGTGGCTCAACCACAGAATGGCCTTCGATACTGGGACGCTTGAGTATGGTAGAGGTTGGTGGAACTGCGAGTGTAGAGGTGAAATTCGTAGATATTCGCAAGAACACCGGTGGCGAAGGCGGCCAACTGGACCATCACTGACGCTGAGGCGCGAAAGCGTGGGGAGCAAACA"                         
##  [9] "TGGGGAATATTGGACAATGGGCGAAAGCCTGATCCAGCAATGCCGCGTGAGTGATGAAGGCCCTAGGGTTGTAAAGCTCTTTTACCCGGGAAGATAATGACTGTACCGGGAGAATAAGCCCCGGCTAACTCCGTGCCAGCAGCCGCGGTAATACGGAGGGGGCTAGCGTTGTTCGGAATTACTGGGCGTAAAGCGCACGTAGGCGGCTTTGTAAGTCAGAGGTGAAAGCCTGGAGCTCAACTCCAGAACTGCCTTTGAGACTGCATCGCTTGAATCCAGGAGAGGTCAGTGGAATTCCGAGTGTAGAGGTGAAATTCGTAGATATTCGGAAGAACACCAGTGGCGAAGGCGGCTGACTGGACTGGTATTGACGCTGAGGTGCGAAAGCGTGGGGAGCAAACA"                         
## [10] "TGGGGAATATTGGACAATGGGCGAAAGCCTGATCCAGCAATGCCGCGTGTGTGAAGAAGGTCTTCGGATTGTAAAGCACTTTCGGCGGGGACGATGATGACGGTACCCGCAGAAGAAGCCCCGGCTAACTTCGTGCCAGCAGCCGCGGTAATACGAAGGGGGCTAGCGTTGCTCGGAATGACTGGGCGTAAAGGGCGCGTAGGCGGTTGCATAAGTTAGATGTGAAATTCCCGGGCTTAACCTGGGGACTGCATTTGATACTATGTGGCTTGAGTATGGAAGAGGGTCGTGGAATTCCCAGTGTAGAGGTGAAATTCGTAGATATTGGGAAGAACACCGGTGGCGAAGGCGGCGACCTGGTCCATAACTGACGCTGAGGCGCGAAAGCGTGGGGAGCAAACA"                         
## [11] "TGGGGAATATTGCACAATGGGCGCAAGCCTGATGCAGCAACGCCGCGTGAGGGACGACGGCCTTCGGGTTGTAAACCTCTTTTAGTAGGGAAGAAGCGAAAGTGACGGTACCTGCAGAAAAAGCACCGGCTAACTACGTGCCAGCAGCCGCGGTAATACGTAGGGTGCAAGCGTTATCCGGAATTATTGGGCGTAAAGAGCTCGTAGGCGGTTTGTCGCGTCTGCTGTGAAAACTGGGGGCTCAACCCCCAGCCTGCAGTGGGTACGGGCAGACTAGAGTGCGGTAGGGGAGATTGGAATTCCTGGTGTAGCGGTGGAATGCGCAGATATCAGGAGGAACACCAATGGCGAAGGCAGATCTCTGGGCCGTAACTGACGCTGAGGAGCGAAAGCATGGGGAGCGAACA"                    
## [12] "TGGGGAATTTTGGACAATGGGCGAAAGCCTGATCCAGCAATGCCGCGTGAGTGAAGAAGGCCTTCGGGTTGTAAAGCTCTTTTGTCAGGGAAGAAACGGTGTGGGCTAATATCCTGCACTAATGACGGTACCTGAAGAATAAGCACCGGCTAACTACGTGCCAGCAGCCGCGGTAATACGTAGGGTGCAAGCGTTAATCGGAATTACTGGGCGTAAAGCGTGCGCAGGCGGTTTTGTAAGTCTGTCGTGAAATCCCCGGGCTCAACCTGGGAATTGCGATGGAGACTGCAAGGCTAGAATCTGGCAGAGGGGGGTAGAATTCCACGTGTAGCAGTGAAATGCGTAGAGATGTGGAGGAACACCGATGGCGAAGGCAGCCCCCTGGGTCAAGATTGACGCTCATGCACGAAAGCGTGGGGAGCAAACA"
## [13] "TGGGGAATATTGGACAATGGGCGAAAGCCTGATCCAGCAATGCCGCGTGAGTGATGAAGGCCTTAGGGTTGTAAAGCTCTTTTACCCGGGATGATAATGACAGTACCGGGAGAATAAGCTCCGGCTAACTCCGTGCCAGCAGCCGCGGTAATACGGAGGGAGCTAGCGTTATTCGGAATTACTGGGCGTAAAGCGCACGTAGGCGGCTTTGTAAGTTAGAGGTGAAAGCCTGGAGCTCAACTCCAGAATTGCCTTTAAGACTGCATCGCTTGAATCCAGGAGAGGTGAGTGGAATTCCGAGTGTAGAGGTGAAATTCGTAGATATTCGGAAGAACACCAGTGGCGAAGGCGGCTCACTGGACTGGTATTGACGCTGAGGTGCGAAAGCGTGGGGAGCAAACA"                         
## [14] "TGGGGAATTTTGGACAATGGGCGAAAGCCTGATCCAGCAATGCCGCGTGCAGGATGAAGGCCTTCGGGTTGTAAACTGCTTTTGTACGGAACGAAAAAGCTCCTTCTAATACAGGGGGCCCATGACGGTACCGTAAGAATAAGCACCGGCTAACTACGTGCCAGCAGCCGCGGTAATACGTAGGGTGCGAGCGTTAATCGGAATTACTGGGCGTAAAGCGTGCGCAGGCGGTTATGTAAGACAGATGTGAAATCCCCGGGCTCAACCTGGGAACTGCATTTGTGACTGCATGGCTAGAGTACGGTAGAGGGGGATGGAATTCCGCGTGTAGCAGTGAAATGCGTAGATATGCGGAGGAACACCGATGGCGAAGGCAATCCCCTGGACCTGTACTGACGCTCATGCACGAAAGCGTGGGGAGCAAACA"
## [15] "TGGGGAATATTGCACAATGGGCGAAAGCCTGATGCAGCGACGCCGCGTGAGGGATGACGGCCTTCGGGTTGTAAACCTCTTTCAGTAGGGAACAAGGCCACACGTGGTGTGGTTGAGGGTACTTGCAGAAGAAGCGCCGGCTAACTACGTGCCAGCAGCCGCGGTAATACGTAGGGCGCAAGCGTTATCCGGAATTATTGGGCGTAAAGAGCTCGTAGGCGGTTTGTCGCGTCTGCCGTGAAAGTCCGGGGCTCAACTCCGGATCTGCGGTGGGTACGGGCAGACTAGAGTGATGTAGGGGAGACTGGAATTCCTGGTGTAGCGGTGAAATGCGCAGATATCAGGAGGAACACCGATGGCGAAGGCAGGTCTCTGGGCATTAACTGACGCTGAGGAGCGAAAGCATGGGGAGCGAACA"         
## [16] "TGGGGAATATTGCACAATGGGCGAAAGCCTGATGCAGCGACGCCGCGTGAGGGATGACGGCCTTCGGGTTGTAAACCTCTTTCAGCAGGGACGAAGCGAGAGTGACGGTACCTGCAGAAGAAGCACCGGCCAACTACGTGCCAGCAGCCGCGGTAATACGTAGGGTGCGAGCGTTGTCCGGAATTACTGGGCGTAAAGAGCTCGTAGGCGGTTTGTCACGTCGGCTGTGAAAACCCATCGCTCAACGGTGGGCTTGCAGTCGATACGGGCTGACTTGAGTACTGCAGGGGAGACTGGAATTCCTGGTGTAGCGGTGAAATGCGCAGATATCAGGAGGAACACCGGTGGCGAAGGCGGGTCTCTGGGCAGTAACTGACGCTGAGGAGCGAAAGCGTGGGGAGCGAACA"                    
## [17] "TGGGGAATATTGCGCAATGGGCGAAAGCCTGACGCAGCGACGCCGCGTGGGGGATGAAGGCCTTCGGGTTGTAAACCTCTTTCAGCTCCGACGAATTCAGACGGTAGGAGCAGAAGAAGCACCGGCCAACTACGTGCCAGCAGCCGCGGTAATACGTAGGGTGCAAGCGTTGTCCGGAATTATTGGGCGTAAAGAGCTCGTAGGCGGTTTGTCGCGTCGGCTGTGAAATCCCGTGGCTCAACTGCGGGTCTGCAGTCGATACGGGCAGACTAGAGTGCGGTAGGGGAGACTGGAATTCCTGGTGTAGCGGTGAAATGCGCAGATATCAGGAGGAACACCGGTGGCGAAGGCGGGTCTCTGGGCCGTAACTGACGCTGAGGAGCGAAAGCGTGGGGAGCGAACA"                        
## [18] "TGGGGAATATTGCACAATGGGCGAAAGCCTGATGCAGCGACGCCGCGTGAGGGATGGAGGCCTTCGGGTTGTAAACCTCTTTCAGCAGGGAAGAAGTGGTCGGGTTCTCTCGGTCGTGACGGTACCTGCAGAAGAAGCGCCGGCTAACTACGTGCCAGCAGCCGCGGTAATACGTAGGGCGCAAGCGTTGTCCGGAATTATTGGGCGTAAAGAGCTCGTAGGCGGTTTGTCGCGTCTGCTGTGAAAACTCGAGGCTCAACCTCGGGCTTGCAGTGGGTACGGGCAGACTAGAGTGCGGTAGGGGAGACTGGAATTCCTGGTGTAGCGGTGGAATGCGCAGATATCAGGAGGAACACCGATGGCGAAGGCAGGTCTCTGGGCCGCAACTGACGCTGAGGAGCGAAAGCATGGGGAGCGAACA"      
## [19] "TGGGGAATATTGCGCAATGGGCGAAAGCCTGACGCAGCGACGCCGCGTGGGGGATGAAGGCCTTCGGGTTGTAAACCTCTTTCAGCTCCGACGAATTCAGACGGTAGGAGCAGAAGAAGCACCGGCCAACTACGTGCCAGCAGCCGCGGTAATACGTAGGGTGCAAGCGTTGTCCGGAATTATTGGGCGTAAAGAGCTCGTAGGCGGTTTGTCGCGTCGGCTGTGAAATCCCGTGGCTCAACTGCGGGTCTGCAGTCGATACGGGCAGACTAGAGTGCGGTAGGGGAGACTGGAATTCCTGGTGTAGCGGTGAAATGCGCAGATATCAGGAGGAACACCGGTGGCGAAGGCGGGTCTCTGGGCCGTAACTGACGCTGAGGAGCGAAAGCGTGGGGAGCAAACA"                        
## [20] "TGGGGAATATTGCGCAATGGGCGAAAGCCTGACGCAGCGACGCCGCGTGAGGGATGACGGCCTTCGGGTTGTAAACCTCTTTCAGCAGGGACGAACACAGACGGTACCTGCAGAAGAAGCACCGGCCAACTACGTGCCAGCAGCCGCGGTAATACGTAGGGTGCAAGCGTTGTCCGGAATTATTGGGCGTAAAGAGCTCGTAGGCGGTCTGTCGCGTCGGCTGTGAAAACTCGGGGCTCAACCCCGAGCCTGCAGTCGATACGGGCAGACTAGAGTGCTGTAGGGGAGACTGGAATTCCTGGTGTAGCGGTGAAATGCGCAGATATCAGGAGGAACACCGGTGGCGAAGGCGGGTCTCTGGGCAGTAACTGACGCTGAGGAGCGAAAGCGTGGGGAGCGAACA"

5a - Genus-level stacked barchart by DNA_or_cDNA~Fox_Aspect (no filter)

5b - Genus-level stacked barchart by DNA_or_cDNA~Snowpack_Layer (filter)

Oh good. But this need improving. Getting the unclassified component of the families for the top main taxa in it, but without a guarantee for them popping into it. Also, basing it on the >0.1% dataset as per Arwyn’s advice.

ps.final.r.glom.genus.psdf_2 <- data.table(speedyseq::psmelt(ps.final.r.glom.genus))
ps.final.r.glom.genus.psdf_2$Genus <- as.character(ps.final.r.glom.genus.psdf_2$Genus)

ps.final.r.glom.genus.psdf_2[, mean := mean(Abundance, na.rm = TRUE), 
                 by = "Genus"]

ps.final.r.glom.genus.psdf_2$Genus[is.na(ps.final.r.glom.genus.psdf_2$Genus) & ps.final.r.glom.genus.psdf_2$Family=="Sphingomonadaceae"] <- "Unc. Sphingomonadaceae"
ps.final.r.glom.genus.psdf_2$Genus[is.na(ps.final.r.glom.genus.psdf_2$Genus) & ps.final.r.glom.genus.psdf_2$Family=="Burkholderiaceae"] <- "Unc. Burkholderiaceae"
ps.final.r.glom.genus.psdf_2$Genus[is.na(ps.final.r.glom.genus.psdf_2$Genus) & ps.final.r.glom.genus.psdf_2$Family=="Labraceae"] <- "Unc. Labraceae"
ps.final.r.glom.genus.psdf_2$Genus[is.na(ps.final.r.glom.genus.psdf_2$Genus) & ps.final.r.glom.genus.psdf_2$Family=="Acetobacteraceae"] <- "Unc. Acetobacteraceae"
ps.final.r.glom.genus.psdf_2$Genus[is.na(ps.final.r.glom.genus.psdf_2$Genus) & ps.final.r.glom.genus.psdf_2$Family=="Microbacteriaceae"] <- "Unc. Microbacteriaceae"
ps.final.r.glom.genus.psdf_2$Genus[is.na(ps.final.r.glom.genus.psdf_2$Genus) & ps.final.r.glom.genus.psdf_2$Family=="Deinococcaceae"] <- "Unc. Deinococcaceae"
ps.final.r.glom.genus.psdf_2$Genus[is.na(ps.final.r.glom.genus.psdf_2$Genus) & ps.final.r.glom.genus.psdf_2$Family=="Gemmatimonaceae"] <- "Unc. Gemmatimonaceae"
ps.final.r.glom.genus.psdf_2$Genus[is.na(ps.final.r.glom.genus.psdf_2$Genus) & ps.final.r.glom.genus.psdf_2$Family=="Intrasporangiaceae"] <- "Unc. Intrasporangiaceae"
ps.final.r.glom.genus.psdf_2$Genus[is.na(ps.final.r.glom.genus.psdf_2$Genus) &          ps.final.r.glom.genus.psdf_2$Family=="Nocardiaceae"] <- "Unc. Nocardiaceae"
ps.final.r.glom.genus.psdf_2$Genus[is.na(ps.final.r.glom.genus.psdf_2$Genus) &          ps.final.r.glom.genus.psdf_2$Family=="Ktedonobacteraceae"] <- "Unc. Ktedonobacteraceae"
ps.final.r.glom.genus.psdf_2$Genus[is.na(ps.final.r.glom.genus.psdf_2$Genus) &          ps.final.r.glom.genus.psdf_2$Family=="Cellulomonadaceae"] <- "Unc. Cellulomonadaceae"
ps.final.r.glom.genus.psdf_2$Genus[is.na(ps.final.r.glom.genus.psdf_2$Genus) &          ps.final.r.glom.genus.psdf_2$Family=="Beijerinckiaceae"] <- "Unc. Beijerinckiaceae"
ps.final.r.glom.genus.psdf_2$Genus[is.na(ps.final.r.glom.genus.psdf_2$Genus) &          is.na(ps.final.r.glom.genus.psdf_2$Family)] <- "Unclassified"
ps.final.r.glom.genus.psdf_2$Genus[is.na(ps.final.r.glom.genus.psdf_2$Genus) &  is.na(ps.final.r.glom.genus.psdf_2$Family) & ps.final.r.glom.genus.psdf_2$Order=="Frankiales"] <- "Unc. Frankiales"

ps.final.r.glom.genus.psdf_2[(mean <= 0.01), Genus := "Other (<1%)"]

ps.final.r.glom.genus.psdf_2$Genus = factor(ps.final.r.glom.genus.psdf_2$Genus,
                         levels = c("Acidiphilium","Arthrobacter",
                                    "Bacillus",
                                    "Bosea",
                                  "Burkholderia-Caballeronia-Paraburkholderia",
                                    "Cryobacterium","Cupriavidus",
                                    "Delftia","Hymenobacter",
                                    "Labrys","Leptothrix",
                                    "Massilia","Methylobacterium",
                                    "Noviherbaspirillum","Oryzihumus",
                                    "Polymorphobacter",
                                    "Salinibacterium",
                                    "Sphingomonas","Other (<1%)"))
tol17rainbow=c("#771155", "#AA4488", "#CC99BB", "#114477",
               "#4477AA", "#77AADD", "#117777", "#44AAAA", 
               "#77CCCC", "#777711", "#AAAA44", "#DDDD77", 
               "#774411", "#AA7744", "#DDAA77", "#771122", 
               "#AA4455", "#DD7788", "azure3")

ggplot(ps.final.r.glom.genus.psdf_2, 
       aes(x = Month, y = Abundance, fill = Genus)) + 
  facet_grid(DNA_or_cDNA~Snowpack_Layer,
             scales = "free",
             space = "free") +
  geom_bar(stat = "identity",
           position = "fill",
           color="black") +
  scale_fill_manual(values = tol17rainbow,
                    na.value = "gray40") +
  scale_y_continuous(expand = c(0,0),
                     labels = scales::percent,
                     breaks = c(0.25,0.5,0.75,1)) +
  theme(axis.title.x = element_blank(),
        axis.text.x = element_text(angle = 35, hjust = 1),
        strip.text.x = element_text(face = "bold",
                                    size = 16),
        strip.background = element_rect(size=1),
        legend.text = element_text(size= 14)) + 
  guides(fill = guide_legend(keywidth = 1, keyheight = 1,
                             nrow = )) +
  ylab("Relative Abundance\n")

5c - Genus-level stacked barchart by DNA_or_cDNA~Fox_Aspect (filter)

ggplot(ps.final.r.glom.genus.psdf_2, 
                                   aes(x = Month, y = Abundance, fill = Genus)) + 
  facet_grid(DNA_or_cDNA~Fox_Aspect,
             scales = "free",
             space = "free") +
  geom_bar(stat = "identity",
           position = "fill",
           color="black") +
  scale_fill_manual(values = tol17rainbow,
                    na.value = "gray40") +
  scale_y_continuous(expand = c(0,0),
                     labels = scales::percent,
                     breaks = c(0.25,0.5,0.75,1)) +
  theme(axis.title.x = element_blank(),
        axis.text.x = element_text(angle = 35, hjust = 1),
        strip.text.x = element_text(face = "bold",
                                    size = 16),
        strip.background = element_rect(size=1),
        legend.text = element_text(size= 14)) + 
  guides(fill = guide_legend(keywidth = 1, keyheight = 1,
                             nrow = )) +
  ylab("Relative Abundance\n")

Lest us not forget that size of Other (<1%) corresponds to the different samples included in that category.

I don’t like this, but plotting all samples now, without grouping.

# ps.final.r.glom.genus.psdf_2$Sample_Name <- as.character(ps.final.r.glom.genus.psdf_2$Sample_Name)
# ps.final.r.glom.genus.psdf_2$Sample_Name <- factor(ps.final.r.glom.genus.psdf_2$Sample_Name,
#                                                  levels=unique(ps.final.r.glom.genus.psdf_2$Sample_Name))
ggplot(ps.final.r.glom.genus.psdf_2,
                                   aes(x = Sample_Name, y = Abundance, fill = Genus)) +
  facet_grid(DNA_or_cDNA~Month,
             scales = "free",
             space = "free") +
  geom_bar(stat = "identity",
           position = "fill",
           color="black") +
  scale_fill_manual(values = tol17rainbow,
                    na.value = "gray40") +
  scale_y_continuous(expand = c(0,0),
                     labels = scales::percent,
                     breaks = c(0.25,0.5,0.75,1)) +
  theme(axis.title.x = element_blank(),
        axis.text.x = element_text(angle = 35, hjust = 1),
        strip.text.x = element_text(face = "bold",
                                    size = 16),
        strip.background = element_rect(size=1),
        legend.text = element_text(size= 14)) +
  guides(fill = guide_legend(keywidth = 1, keyheight = 1,
                             nrow = )) +
  ylab("Relative Abundance\n")

6 - Unconstrained ordination

nmds_shapes = c(15,16,17,18,6)
nmdscolors_by_month = c("April" = "#A6048B",
                        "June" = "#0704A6",
                        "Early July" = "#727A0B",
                        "Late July" = "#6C0504")
nmdscolors_by_snpl = c("Top" = "#151AC6",
                       "Mid" = "#C6156D",
                       "Basal" = "#890909",
                       "Sup. Ice" = "#091D89",
                       "Glacial Ice" = "#580303")

ps.final.r.dna = subset_samples(ps.final.r, DNA_or_cDNA != "cDNA")
ps.final.r.cdna = subset_samples(ps.final.r, DNA_or_cDNA != "DNA")

ps.final.r.dna.nmds  <- ordinate(ps.final.r.dna,
                         "NMDS",
                         distance="jsd",
                         maxit = 1e3)
## Run 0 stress 0.2107623 
## Run 1 stress 0.2270207 
## Run 2 stress 0.2259277 
## Run 3 stress 0.2158645 
## Run 4 stress 0.2249122 
## Run 5 stress 0.2276945 
## Run 6 stress 0.2254988 
## Run 7 stress 0.2087417 
## ... New best solution
## ... Procrustes: rmse 0.1326501  max resid 0.2782618 
## Run 8 stress 0.2034057 
## ... New best solution
## ... Procrustes: rmse 0.0775235  max resid 0.4113966 
## Run 9 stress 0.2030457 
## ... New best solution
## ... Procrustes: rmse 0.1023723  max resid 0.4191121 
## Run 10 stress 0.2178268 
## Run 11 stress 0.2064812 
## Run 12 stress 0.2106171 
## Run 13 stress 0.2087421 
## Run 14 stress 0.2030472 
## ... Procrustes: rmse 0.0006936713  max resid 0.00327211 
## ... Similar to previous best
## Run 15 stress 0.2264449 
## Run 16 stress 0.2019616 
## ... New best solution
## ... Procrustes: rmse 0.04318994  max resid 0.188176 
## Run 17 stress 0.208532 
## Run 18 stress 0.2236278 
## Run 19 stress 0.2161764 
## Run 20 stress 0.2188811 
## *** No convergence -- monoMDS stopping criteria:
##      1: no. of iterations >= maxit
##     18: stress ratio > sratmax
##      1: scale factor of the gradient < sfgrmin
ps.final.r.cdna.nmds  <- ordinate(ps.final.r.cdna,
                         "NMDS",
                         distance="jsd",
                         maxit = 1e3)
## Run 0 stress 0.1654426 
## Run 1 stress 0.1595584 
## ... New best solution
## ... Procrustes: rmse 0.1139479  max resid 0.3178704 
## Run 2 stress 0.1710025 
## Run 3 stress 0.1718853 
## Run 4 stress 0.1648546 
## Run 5 stress 0.1678885 
## Run 6 stress 0.1649862 
## Run 7 stress 0.1686493 
## Run 8 stress 0.1656653 
## Run 9 stress 0.166164 
## Run 10 stress 0.1599933 
## ... Procrustes: rmse 0.0790213  max resid 0.2159161 
## Run 11 stress 0.162065 
## Run 12 stress 0.1724002 
## Run 13 stress 0.176513 
## Run 14 stress 0.1667006 
## Run 15 stress 0.1616702 
## Run 16 stress 0.1697515 
## Run 17 stress 0.1661504 
## Run 18 stress 0.1651428 
## Run 19 stress 0.1605854 
## Run 20 stress 0.1687601 
## *** No convergence -- monoMDS stopping criteria:
##     19: stress ratio > sratmax
##      1: scale factor of the gradient < sfgrmin
ps.final.r.dna.nmds.plot <- plot_ordination(ps.final.r.dna, ps.final.r.dna.nmds) +
               geom_point(aes(colour=Month, shape=Snowpack_Layer),
                          size=5, alpha = 0.5, stroke=1.5) +
               scale_shape_manual(values=nmds_shapes) +
               scale_colour_manual(values=nmdscolors_by_month) +
               scale_fill_manual(values=nmdscolors_by_month) +
               stat_ellipse(aes(fill=Month),
                          geom = "polygon",
                          type="t",
                          alpha=0.15,
                          linetype = 2,
                          show.legend = FALSE) +
  ggtitle("DNA nMDS") +
  guides(colour = guide_legend(title="Month"))

ps.final.r.cdna.nmds.plot <- plot_ordination(ps.final.r.cdna, ps.final.r.cdna.nmds) +
               geom_point(aes(colour=Month, shape=Snowpack_Layer),
                          size=5, alpha = 0.7, stroke=1.5) +
               scale_shape_manual(values=nmds_shapes) +
               scale_colour_manual(values=nmdscolors_by_month) +
               scale_fill_manual(values=nmdscolors_by_month) +
               stat_ellipse(aes(fill=Month),
                          geom = "polygon",
                          type="t",
                          alpha=0.2,
                          linetype = 2,
                          show.legend = FALSE)+
  ggtitle("cDNA nMDS") +
  guides(colour = guide_legend(title="Month"))
plot_grid(ps.final.r.dna.nmds.plot + theme_bw(), 
          ps.final.r.cdna.nmds.plot + theme_bw())

7 - Alpha-diversity metrics.

month_comparisons = list( c("April", "June"), c("April", "Early July"), 
                          c("April", "Late July"),
                          c("June", "Early July"), c("June", "Late July"),
                          c("Early July", "Late July"))

site_comparisons = list( c("NW_AWS", "N_NE"),
                         c("NW_AWS", "SWS1SE"),
                         c("N_NE", "SWS1SE"))

slayer_comparisons = list( c("TOP", "MID"), c("TOP", "BASAL"), c("TOP", "SUP ICE"), 
                           c("TOP", "GL ICE"),
                           c("MID", "BASAL"), c("MID", "SUP ICE"), 
                           c("MID", "GL ICE"),
                           c("BASAL", "SUP ICE"), c("BASAL", "GL ICE"),
                           c("SUP ICE", "GL ICE"))

All measures for Sample_ID, no stats, facet by DNA/cDNA

plot_richness(ps.final,
              x="Sample_ID", color="Sample_ID",
              measures=c("Observed","Shannon","Simpson")) +
  geom_jitter(shape=16, position=position_jitter(0.2)) + 
  facet_wrap(DNA_or_cDNA~variable, scales = "free") +
  labs(color="Sample_ID") +
  guides(color=guide_legend(ncol=2)) +
  theme(axis.text.x = element_text(angle = 20,
                                   vjust = 1,
                                   hjust = 1,
                                   size = 12, face = "bold"),
        axis.text.y = element_text(face = "bold"),
        axis.title = element_blank(),
        strip.text = element_text(size = 12, margin = margin()))

All measures by Sample_ID

mini_metadata = metadata %>% 
  select(c(Sample_ID, Sample_Name, Month, DNA_or_cDNA, Fox_Aspect, Snowpack_Layer))

ps.final.measures = estimate_richness(ps.final, measures = c("Observed","Shannon","Simpson")) %>% 
  rownames_to_column(var = "Sample_ID") %>% 
  mutate(Sample_ID = gsub("\\.", "-", Sample_ID))

full_join(ps.final.measures, mini_metadata, 
                                       by="Sample_ID") %>% 
  kable(digits = 4) %>%
  kable_styling(full_width = T, bootstrap_options = c("striped", "hover")) %>% 
  row_spec(0, bold=TRUE)
Sample_ID Observed Shannon Simpson Sample_Name Month DNA_or_cDNA Fox_Aspect Snowpack_Layer
A03-cDNA-1 286 0.7114 0.2103 April_N_NE_Basal_R April cDNA N_NE BASAL
A03-DNA-1 138 3.2774 0.9407 June_N_NE_Top20_D June DNA N_NE TOP
A04-cDNA-2 163 2.0361 0.7852 April_NW_AWS_Top_R April cDNA NW_AWS TOP
A04-DNA-2 75 1.6409 0.6083 April_ NW_AWS_ Basal_D April DNA NW_AWS BASAL
A05-cDNA-3 515 2.2067 0.7048 June_NW_AWS_Basal_R June cDNA NW_AWS BASAL
A05-DNA-3 1390 5.1230 0.9813 30 July_ N_NE_ Gl ice_D Late July DNA N_NE GL ICE
A06-cDNA-4 125 1.3510 0.6058 June_N_NE_Basal_R June cDNA N_NE BASAL
A06-DNA-4 210 3.9337 0.9650 June_ SWS1SE_ Top20_D June DNA SWS1SE TOP
A07-cDNA-5 1635 4.6591 0.9488 9 July_NW_AWS_Mid_R Early July cDNA NW_AWS MID
A07-DNA-5 609 3.3602 0.8726 9 July_ NW_AWS_ Sup ice_D Early July DNA NW_AWS SUP ICE
A08-DNA-6 563 4.0141 0.8786 April_ NW_AWS_ Mid_D April DNA NW_AWS MID
A09-cDNA-7 1875 5.5902 0.9807 30 July_NW_AWS_Top20_R Late July cDNA NW_AWS TOP
A10-cDNA-8 842 4.7339 0.9520 30 July_N_NE_Gl ice_R Late July cDNA N_NE GL ICE
A10-DNA-8 179 3.3916 0.9493 April_NW_AWS_Top20_D April DNA NW_AWS TOP
A11-cDNA-9 3 0.2518 0.1121 April_NW_AWS_Mid_R April cDNA NW_AWS MID
A11-DNA-9 1819 5.3121 0.9764 9 July_SWS1SE_Gl ice 1_D Early July DNA SWS1SE GL ICE
A12-cDNA-10 54 1.4929 0.6133 June_NW_AWS_Top20_R June cDNA NW_AWS TOP
A12-DNA-10 2719 5.7179 0.9837 30 July_N_NE_Top20_D Late July DNA N_NE TOP
B01-cDNA-11 453 4.0612 0.9622 30 July_NW_AWS_Sup ice_R Late July cDNA NW_AWS SUP ICE
B01-DNA-11 187 3.5277 0.9546 June_NW_AWS_Mid_D June DNA NW_AWS MID
B02-cDNA-12 733 4.5851 0.9683 30 July_N_NE_Top20_R Late July cDNA N_NE TOP
B02-DNA-12 310 2.1151 0.7154 9 July_N_NE_Mid_D Early July DNA N_NE MID
B03-cDNA-13 38 1.1241 0.4091 June_SWS1SE_Basal_R June cDNA SWS1SE BASAL
B03-DNA-13 457 4.2859 0.9559 April_SWS1SE_Mid_D April DNA SWS1SE MID
B04-cDNA-14 67 1.1867 0.4139 9 July_NW_AWS_Top20_R Early July cDNA NW_AWS TOP
B05-cDNA-15 1307 5.7881 0.9881 9 July_NW_AWS_Gl ice 1_R Early July cDNA NW_AWS GL ICE
B06-cDNA-16 134 1.4982 0.6595 April_SWS1SE_Top20_R April cDNA SWS1SE TOP
B06-DNA-2 156 0.6586 0.2271 April_SWS1SE_Basal_D April DNA SWS1SE BASAL
B07-DNA-3 1302 5.3971 0.9910 June_NW AWS_Top20_D June DNA NW_AWS TOP
B08-cDNA-18 21 0.9250 0.5154 June_NW_AWS_Mid_R June cDNA NW_AWS MID
B08-DNA-4 986 5.2137 0.9805 30 July_SWS1SE_Gl ice_D Late July DNA SWS1SE GL ICE
B09-DNA-5 105 2.9238 0.9150 June_N_NE_Mid_D June DNA N_NE MID
B10-cDNA-21 400 4.5700 0.9771 9 July_SWS1SE_Sup ice_R Early July cDNA SWS1SE SUP ICE
B10-DNA-6 377 4.1547 0.9598 June_SWS1SE_Basal_D June DNA SWS1SE BASAL
B12-cDNA-23 29 0.9625 0.3511 30 July_SWS1SE_Slush_R Late July cDNA SWS1SE TOP
B12-DNA-8 1998 5.3382 0.9817 9 July_ NW_AWS_Gl ice 1_D Early July DNA NW_AWS GL ICE
C01-DNA-9 165 3.3370 0.9306 June_SWS1SE_Mid_D June DNA SWS1SE MID
C02-cDNA-26 907 4.1753 0.9419 30 July_SWS1SE_Gl ice_R Late July cDNA SWS1SE GL ICE
C03-cDNA-27 537 3.3393 0.8840 9 July_N_NE_Sup ice_R Early July cDNA N_NE SUP ICE
C03-DNA-B2-1 677 5.1429 0.9887 June_NW_AWS_Basal_D June DNA NW_AWS BASAL
C04-cDNA-28 307 3.3747 0.8749 30 July_NW_AWS_Gl ice_R Late July cDNA NW_AWS GL ICE
C04-DNA-2 715 4.2502 0.9430 30 July_NW_AWS_Gl ice_D Late July DNA NW_AWS GL ICE
C05-DNA-3 1739 5.2347 0.9838 30 July_NW_AWS_Sup ice_D Late July DNA NW_AWS SUP ICE
C06-DNA-4 236 3.4402 0.9006 April_SWS1SE_Top20_D April DNA SWS1SE TOP
C07-DNA-5 2789 5.9749 0.9894 9 July_NW_AWS_Mid_D Early July DNA NW_AWS MID
C08-DNA-6 307 4.5001 0.9837 June_N_NE_Basal_D June DNA N_NE BASAL
C09-cDNA-33 36 1.3284 0.5802 June_N_NE_Top20_R June cDNA N_NE TOP
C09-DNA-7 148 3.8243 0.9690 9 July_N_NE_Gl ice_D Early July DNA N_NE GL ICE
C10-DNA-8 383 3.2676 0.8801 9 July_ SWS1SE_Mid_D Early July DNA SWS1SE MID
C11-cDNA-35 65 1.5747 0.6657 9 July_SWS1SE_Mid_R Early July cDNA SWS1SE MID
C12-cDNA-36 19 0.4213 0.1462 9 July_N_NE_Gl ice_R Early July cDNA N_NE GL ICE
C12-DNA-10 396 2.2692 0.7492 9 July_N_NE_Sup ice_D Early July DNA N_NE SUP ICE
D01-DNA-B4-1 147 0.6873 0.2422 30 July_SWS1SE_Slush_D Late July DNA SWS1SE TOP
D02-cDNA-38 128 2.1030 0.7780 April_SWS1SE_Mid_R April cDNA SWS1SE MID
D03-cDNA-39 266 3.4007 0.9035 April_N_NE_Mid_R April cDNA N_NE MID
D04-cDNA-40 79 1.5083 0.6114 June_SWS1SE_Top20_R June cDNA SWS1SE TOP
D05-cDNA-41 154 1.8281 0.7377 June_N_NE_Mid_R June cDNA N_NE MID
D05-DNA-5 90 0.6634 0.2169 9 July_N_NE_Top20_D Early July DNA N_NE TOP
D06-cDNA-42 841 3.6042 0.9024 9 July_NW_AWS_Sup ice_R Early July cDNA NW_AWS SUP ICE
D06-DNA-6 2292 5.5649 0.9810 9 July_SWS1SE_Sup ice_D Early July DNA SWS1SE SUP ICE
D07-cDNA-43 1444 4.5026 0.9489 9 July_SWS1SE_Gl ice 2_R Early July cDNA SWS1SE GL ICE
D07-DNA-7 467 1.8746 0.5308 9 July_NW_AWS_Top20_D Early July DNA NW_AWS TOP
D08-cDNA-44 120 2.6011 0.8428 9 July_N_NE_Mid_R Early July cDNA N_NE MID
D08-DNA-8 1045 4.1510 0.9245 30 July_NW_AWS_Top20_D Late July DNA NW_AWS TOP
D09-DNA-9 2986 5.9104 0.9895 9 July_SWS1SE_Top20_D Early July DNA SWS1SE TOP
B05-DNA-B1-1 NA NA NA April_N_NE_Top20_D April DNA N_NE_D TOP
C02-DNA-10 NA NA NA April_N_NE_Mid_D April DNA N_NE_D MID
C11-DNA-9 NA NA NA April_N_NE_Basal_D April DNA N_NE_D BASAL
D02-DNA-2 NA NA NA 9 July_NW_AWS_Gl ice 2_D Early July DNA NW_AWS GL ICE
B11-DNA-7 NA NA NA 9 July_SWS1SE_Gl ice 2_D Early July DNA SWS1SE GL ICE
A01-DNA-MM NA NA NA MM_D PCR_Control DNA Control Control
A02-DNA-MM-H20 NA NA NA MM+H2O_D PCR_Control DNA Control Control
A09-DNA-7 NA NA NA St MQ_D Control DNA Control Control
B04-DNA-14 NA NA NA St MQ Unis_D Control DNA Control Control
D03-DNA-3 NA NA NA SLS_D Control DNA Control Control
D04-DNA-4 NA NA NA ST Control_D Control DNA Control Control
D10-DNA-B5-1 NA NA NA St MQ_D Control DNA Control Control
D11-DNA-2 NA NA NA SLS St_D Control DNA Control Control
D12-DNA-3 NA NA NA St Only_D Control DNA Control Control
E01-DNA-4 NA NA NA St MQ2_D Control DNA Control Control
E02-DNA-MM NA NA NA MM_D PCR_Control DNA Control Control
E03-DNA-MM-H20 NA NA NA MM+H20_D PCR_Control DNA Control Control
C07-cDNA-31 NA NA NA April_NW_AWS_Basal_R April cDNA NW_AWS BASAL
B07-cDNA-17 NA NA NA April_N_NE_Top20_R April cDNA N_NE TOP
C08-cDNA-32 NA NA NA April_SWS1SE_Basal_R April cDNA SWS1SE BASAL
B09-cDNA-19 NA NA NA June_SWS1SE_Mid_R June cDNA SWS1SE MID
C10-cDNA-34 NA NA NA 9 July_NW_AWS_Gl ice 2_R Early July cDNA NW_AWS GL ICE
B11-cDNA-22 NA NA NA 9 July_ N_NE_Top20_R Early July cDNA N_NE TOP
A08-cDNA-6 NA NA NA 9 July_SWS1SE_Gl ice_R Early July cDNA SWS1SE GL ICE
A01-cDNA-MM NA NA NA MM_R PCR_Control cDNA Control Control
A02-cDNA-MM-H20 NA NA NA MM+H20_R PCR_Control cDNA Control Control
C01-cDNA-25 NA NA NA St MQ Unis_R Control cDNA Control Control
C05-cDNA-29 NA NA NA SLS St_R Control cDNA Control Control
C06-cDNA-30 NA NA NA St MQ2_R Control cDNA Control Control
D01-cDNA-37 NA NA NA St Only_R Control cDNA Control Control
D09-cDNA-45 NA NA NA SLS St_R Control cDNA Control Control
D10-cDNA-46 NA NA NA St MQ_R Control cDNA Control Control
D11-cDNA-47 NA NA NA St Only_R Control cDNA Control Control
D12-cDNA-MM NA NA NA MM_R PCR_Control cDNA Control Control
E01-cDNA-MM-H20 NA NA NA MM+H20_R PCR_Control cDNA Control Control

All measures for Month, stats included, facet by DNA/cDNA

plot_richness(ps.final,
              x="Month", color="Month",
              measures=c("Observed","Shannon","Simpson")) +
  geom_violin(alpha = .5) +
  geom_jitter(shape=16, position=position_jitter(0.2)) + 
  facet_wrap(DNA_or_cDNA~variable, scales = "free") +
  # stat_compare_means(comparisons = my_comparisons)+
  stat_compare_means(comparisons = month_comparisons,
                     method = "t.test",
                     var.equal="welch",
                     label.x.npc = "right", 
                     label.y.npc = "top",
                     symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1), 
                     symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
  labs(color="Month") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 20,
                                   vjust = 1,
                                   hjust = 1,
                                   size = 12, face = "bold"),
        axis.text.y = element_text(face = "bold"),
        axis.title = element_blank(),
        strip.text = element_text(size = 12, margin = margin()))

All measures for Month, stats included, facet by DNA/cDNA, “GL ICE” removed, “Gl ice” samples removed, samples “30 July_SWS1SE_Slush_R”, “30 July_SWS1SE_Slush_D” also removed

plot_richness(subset_samples(ps.final, 
                             Snowpack_Layer != "GL ICE" & 
                             Sample_Name != "30 July_SWS1SE_Slush_R" &
                             Sample_Name != "30 July_SWS1SE_Slush_D" &
                             !(Sample_Name %in% gl_ice_spls)),
              x="Month", color="Month",
              measures=c("Observed","Shannon","Simpson")) +
  geom_violin(alpha = .5) +
  geom_jitter(shape=16, position=position_jitter(0.2)) + 
  facet_wrap(DNA_or_cDNA~variable, scales = "free") +
  # stat_compare_means(comparisons = my_comparisons)+
  stat_compare_means(comparisons = month_comparisons,
                     method = "t.test",
                     var.equal="welch",
                     label.x.npc = "right", 
                     label.y.npc = "top",
                     symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1), 
                     symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
  labs(color="Month") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 20,
                                   vjust = 1,
                                   hjust = 1,
                                   size = 12, face = "bold"),
        axis.text.y = element_text(face = "bold"),
        axis.title = element_blank(),
        strip.text = element_text(size = 12, margin = margin()))

All measures for Month, stats included, facet by DNA/cDNA, no (“Early July”, “Late July”)

mini_month_comparisons = list(c("April", "June"), 
                          c("April", "June"))

plot_richness(subset_samples(ps.final, Month != "Early July" & Month != "Late July"),
              x="Month", color="Month",
              measures=c("Observed","Shannon","Simpson")) +
  geom_violin(alpha = .5) +
  geom_jitter(shape=16, position=position_jitter(0.2)) + 
  facet_wrap(DNA_or_cDNA~variable, scales = "free") +
  # stat_compare_means(comparisons = my_comparisons)+
  stat_compare_means(comparisons = mini_month_comparisons,
                     method = "t.test",
                     var.equal="welch",
                     label.x.npc = "right", 
                     label.y.npc = "top",
                     symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1), 
                     symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
  labs(color="Month") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 20,
                                   vjust = 1,
                                   hjust = 1,
                                   size = 12, face = "bold"),
        axis.text.y = element_text(face = "bold"),
        axis.title = element_blank(),
        strip.text = element_text(size = 12, margin = margin()))

All measures for Fox_Aspect, no stats, facet by DNA/cDNA

#now without stats, only shapes
plot_richness(ps.final,
              x="Fox_Aspect", color="Fox_Aspect",
              measures=c("Observed","Shannon","Simpson")) +
  geom_violin() +
  geom_jitter(aes(shape=Month), size=5, alpha=0.7,
              position=position_jitter(0.2)) + 
  facet_wrap(DNA_or_cDNA~variable, scales = "free") +
  labs(color="Fox Aspect") +
  scale_shape_manual(values=nmds_shapes) +
  theme(axis.text.x = element_text(angle = 20,
                                   vjust = 1,
                                   hjust = 1,
                                   size = 12, face = "bold"),
        axis.text.y = element_text(face = "bold"),
        axis.title = element_blank(),
        strip.text = element_text(size = 12, margin = margin()))

All measures for Fox_Aspect, stats included, facet by DNA/cDNA

plot_richness(ps.final,
              x="Fox_Aspect", color="Fox_Aspect",
              measures=c("Observed","Shannon","Simpson")) +
  geom_violin(alpha=0.5) +
  geom_jitter(shape=16, position=position_jitter(0.2)) + 
  facet_wrap(DNA_or_cDNA~variable, scales = "free") +
  # stat_compare_means(comparisons = my_comparisons)+
  stat_compare_means(comparisons = site_comparisons,
                     method = "t.test",
                     var.equal="welch",
                     label.x.npc = "right", 
                     label.y.npc = "top",
                     symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1), 
                     symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
  labs(color="Fox Aspect") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 20,
                                   vjust = 1,
                                   hjust = 1,
                                   size = 12, face = "bold"),
        axis.text.y = element_text(face = "bold"),
        axis.title = element_blank(),
        strip.text = element_text(size = 12, margin = margin()))

All measures for Fox_Aspect, stats included, facet by DNA/cDNA, “GL ICE” removed, “Gl ice” samples removed, samples “30 July_SWS1SE_Slush_R”, “30 July_SWS1SE_Slush_D” also removed

plot_richness(subset_samples(ps.final, 
                             Snowpack_Layer != "GL ICE" & 
                             Sample_Name != "30 July_SWS1SE_Slush_R" &
                             Sample_Name != "30 July_SWS1SE_Slush_D" &
                             !(Sample_Name %in% gl_ice_spls)),
              x="Fox_Aspect", color="Fox_Aspect",
              measures=c("Observed","Shannon","Simpson")) +
  geom_violin(alpha=0.5) +
  geom_jitter(shape=16, position=position_jitter(0.2)) + 
  facet_wrap(DNA_or_cDNA~variable, scales = "free") +
  # stat_compare_means(comparisons = my_comparisons)+
  stat_compare_means(comparisons = site_comparisons,
                     method = "t.test",
                     var.equal="welch",
                     label.x.npc = "right", 
                     label.y.npc = "top",
                     symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1), 
                     symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
  labs(color="Fox Aspect") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 20,
                                   vjust = 1,
                                   hjust = 1,
                                   size = 12, face = "bold"),
        axis.text.y = element_text(face = "bold"),
        axis.title = element_blank(),
        strip.text = element_text(size = 12, margin = margin()))

All measures for Fox_Aspect, stats included, facet by DNA/cDNA, no (“Early July”, “Late July”)

plot_richness(subset_samples(ps.final, Snowpack_Layer != "GL ICE" &
                             Month != "Early July" & Month != "Late July"),
              x="Fox_Aspect", color="Fox_Aspect",
              measures=c("Observed","Shannon","Simpson")) +
  geom_violin(alpha = .5) +
  geom_jitter(shape=16, position=position_jitter(0.2)) + 
  facet_wrap(DNA_or_cDNA~variable, scales = "free") +
  # stat_compare_means(comparisons = my_comparisons)+
  stat_compare_means(comparisons = mini_month_comparisons,
                     method = "t.test",
                     var.equal="welch",
                     label.x.npc = "right", 
                     label.y.npc = "top",
                     symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1), 
                     symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
  labs(color="Fox_Aspect") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 20,
                                   vjust = 1,
                                   hjust = 1,
                                   size = 12, face = "bold"),
        axis.text.y = element_text(face = "bold"),
        axis.title = element_blank(),
        strip.text = element_text(size = 12, margin = margin()))
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed

## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed

## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed

## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed

## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed

## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed

All measures for Fox_Aspect, stats included, facet by DNA/cDNA, July only

plot_richness(subset_samples(ps.final, Month == c("Early July", "Late July")),
              x="Fox_Aspect", color="Fox_Aspect",
              measures=c("Observed","Shannon","Simpson")) +
  geom_violin(alpha = .5) +
  geom_jitter(shape=16, position=position_jitter(0.2)) + 
  facet_wrap(DNA_or_cDNA~variable, scales = "free") +
  # stat_compare_means(comparisons = my_comparisons)+
  stat_compare_means(comparisons = mini_month_comparisons,
                     method = "t.test",
                     var.equal="welch",
                     label.x.npc = "right", 
                     label.y.npc = "top",
                     symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1), 
                     symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
  labs(color="Fox_Aspect") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 20,
                                   vjust = 1,
                                   hjust = 1,
                                   size = 12, face = "bold"),
        axis.text.y = element_text(face = "bold"),
        axis.title = element_blank(),
        strip.text = element_text(size = 12, margin = margin()))
## Warning in `==.default`(Month, c("Early July", "Late July")): longer object
## length is not a multiple of shorter object length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed

## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed

## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed

## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed

## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed

## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed

All measures for Fox_Aspect, stats included, facet by DNA/cDNA, July only, “GL ICE” removed, “Gl ice” samples removed, samples “30 July_SWS1SE_Slush_R”, “30 July_SWS1SE_Slush_D” also removed."

plot_richness(subset_samples(ps.final, Month == c("Early July", "Late July"), 
                             Snowpack_Layer != "GL ICE" & 
                             Sample_Name != "30 July_SWS1SE_Slush_R" &
                             Sample_Name != "30 July_SWS1SE_Slush_D" &
                             !(Sample_Name %in% gl_ice_spls)),
              x="Fox_Aspect", color="Fox_Aspect",
              measures=c("Observed","Shannon","Simpson")) +
  geom_violin(alpha = .5) +
  geom_jitter(shape=16, position=position_jitter(0.2)) + 
  facet_wrap(DNA_or_cDNA~variable, scales = "free") +
  # stat_compare_means(comparisons = my_comparisons)+
  stat_compare_means(comparisons = mini_month_comparisons,
                     method = "t.test",
                     var.equal="welch",
                     label.x.npc = "right", 
                     label.y.npc = "top",
                     symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1), 
                     symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
  labs(color="Fox_Aspect") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 20,
                                   vjust = 1,
                                   hjust = 1,
                                   size = 12, face = "bold"),
        axis.text.y = element_text(face = "bold"),
        axis.title = element_blank(),
        strip.text = element_text(size = 12, margin = margin()))
## Warning in `==.default`(Month, c("Early July", "Late July")): longer object
## length is not a multiple of shorter object length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed

## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed

## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed

## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed

## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed

## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed

All measures for Fox_Aspect, stats included, facet by DNA/cDNA, “Early July” only, “GL ICE” removed, “Gl ice” samples removed

plot_richness(subset_samples(ps.final, Snowpack_Layer != "GL ICE" &
                               Month == "Early July" & !(Sample_Name %in% gl_ice_spls)),
              x="Fox_Aspect", color="Fox_Aspect",
              measures=c("Observed","Shannon","Simpson")) +
  geom_violin(alpha = .5) +
  geom_jitter(shape=16, position=position_jitter(0.2)) + 
  facet_wrap(DNA_or_cDNA~variable, scales = "free") +
  # stat_compare_means(comparisons = my_comparisons)+
  stat_compare_means(comparisons = mini_month_comparisons,
                     method = "t.test",
                     var.equal="welch",
                     label.x.npc = "right", 
                     label.y.npc = "top",
                     symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1), 
                     symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
  labs(color="Fox_Aspect") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 20,
                                   vjust = 1,
                                   hjust = 1,
                                   size = 12, face = "bold"),
        axis.text.y = element_text(face = "bold"),
        axis.title = element_blank(),
        strip.text = element_text(size = 12, margin = margin()))
## Warning in max(data$density): no non-missing arguments to max; returning -Inf
## Warning: Computation failed in `stat_ydensity()`:
## replacement has 1 row, data has 0
## Warning in max(data$density): no non-missing arguments to max; returning -Inf
## Warning: Computation failed in `stat_ydensity()`:
## replacement has 1 row, data has 0
## Warning in max(data$density): no non-missing arguments to max; returning -Inf
## Warning: Computation failed in `stat_ydensity()`:
## replacement has 1 row, data has 0
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed

## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed

## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed

## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed

## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed

## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed

All measures for Fox_Aspect, stats included, facet by DNA/cDNA, “Late July” only, “GL ICE”, samples “30 July_SWS1SE_Slush_R”, “30 July_SWS1SE_Slush_D” removed, “Gl ice” samples removed

plot_richness(subset_samples(ps.final, Snowpack_Layer != "GL ICE" & Month == "Late July" &
                               Sample_Name != "30 July_SWS1SE_Slush_R" &
                               Sample_Name != "30 July_SWS1SE_Slush_D" & 
                               !(Sample_Name %in% gl_ice_spls)),
              x="Fox_Aspect", color="Fox_Aspect",
              measures=c("Observed","Shannon","Simpson")) +
  geom_violin(alpha = .5) +
  geom_jitter(shape=16, position=position_jitter(0.2)) + 
  facet_wrap(DNA_or_cDNA~variable, scales = "free") +
  # stat_compare_means(comparisons = my_comparisons)+
  stat_compare_means(comparisons = mini_month_comparisons,
                     method = "t.test",
                     var.equal="welch",
                     label.x.npc = "right", 
                     label.y.npc = "top",
                     symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1), 
                     symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
  labs(color="Fox_Aspect") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 20,
                                   vjust = 1,
                                   hjust = 1,
                                   size = 12, face = "bold"),
        axis.text.y = element_text(face = "bold"),
        axis.title = element_blank(),
        strip.text = element_text(size = 12, margin = margin()))
## Warning in max(data$density): no non-missing arguments to max; returning -Inf
## Warning: Computation failed in `stat_ydensity()`:
## replacement has 1 row, data has 0
## Warning in max(data$density): no non-missing arguments to max; returning -Inf
## Warning: Computation failed in `stat_ydensity()`:
## replacement has 1 row, data has 0
## Warning in max(data$density): no non-missing arguments to max; returning -Inf
## Warning: Computation failed in `stat_ydensity()`:
## replacement has 1 row, data has 0
## Warning in max(data$density): no non-missing arguments to max; returning -Inf
## Warning: Computation failed in `stat_ydensity()`:
## replacement has 1 row, data has 0
## Warning in max(data$density): no non-missing arguments to max; returning -Inf
## Warning: Computation failed in `stat_ydensity()`:
## replacement has 1 row, data has 0
## Warning in max(data$density): no non-missing arguments to max; returning -Inf
## Warning: Computation failed in `stat_ydensity()`:
## replacement has 1 row, data has 0
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed

## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed

## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed

## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed

## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed

## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed

All measures for Snowpack_Layer, no stats, facet by DNA/cDNA

#now without stats, only shapes
plot_richness(ps.final,
              x="Snowpack_Layer", color="Snowpack_Layer",
              measures=c("Observed","Shannon","Simpson")) +
  geom_violin() +
  geom_jitter(aes(shape=Month), size=5, alpha=0.7,
              position=position_jitter(0.2)) + 
  facet_wrap(DNA_or_cDNA~variable, scales = "free") +
  labs(color="Snowpack Layer") +
  scale_shape_manual(values=nmds_shapes) +
  theme(axis.text.x = element_text(angle = 20,
                                   vjust = 1,
                                   hjust = 1,
                                   size = 12, face = "bold"),
        axis.text.y = element_text(face = "bold"),
        axis.title = element_blank(),
        strip.text = element_text(size = 12, margin = margin()))

All measures for Snowpack_Layer, stats included, facet by DNA/cDNA

plot_richness(ps.final,
              x="Snowpack_Layer", color="Snowpack_Layer",
              measures=c("Observed","Shannon","Simpson")) +
  geom_violin(alpha = .5) +
  geom_jitter(position=position_jitter(0.2)) + 
  facet_wrap(DNA_or_cDNA~variable, scales = "free") +
  # stat_compare_means(comparisons = my_comparisons)+
  stat_compare_means(comparisons = slayer_comparisons,
                     method = "t.test",
                     var.equal="welch",
                     label.x.npc = "right", 
                     label.y.npc = "top",
                     symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1), 
                     symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
  labs(color="Snowpack Layer") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 20,
                                   vjust = 1,
                                   hjust = 1,
                                   size = 12, face = "bold"),
        axis.text.y = element_text(face = "bold"),
        axis.title = element_blank(),
        strip.text = element_text(size = 12, margin = margin()))

All measures for Snowpack_Layer in July only (“Early July”, “Late July”), stats included, facet by DNA/cDNA

plot_richness(subset_samples(ps.final, Month = c("Early July", "Late July")),
              x="Snowpack_Layer", color="Snowpack_Layer",
              measures=c("Observed","Shannon","Simpson")) +
  geom_violin(alpha = .5) +
  geom_jitter(position=position_jitter(0.2)) + 
  facet_wrap(DNA_or_cDNA~variable, scales = "free") +
  # stat_compare_means(comparisons = my_comparisons)+
  stat_compare_means(comparisons = slayer_comparisons,
                     method = "t.test",
                     var.equal="welch",
                     label.x.npc = "right", 
                     label.y.npc = "top",
                     symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1), 
                       symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
  labs(color="Snowpack Layer") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 20,
                                   vjust = 1,
                                   hjust = 1,
                                   size = 12, face = "bold"),
        axis.text.y = element_text(face = "bold"),
        axis.title = element_blank(),
        strip.text = element_text(size = 12, margin = margin()))

# ggsave("july_alpha.svg", width = 12, height = 8)
# ggsave("july_alpha.tiff", width = 12, height = 8)

All measures for Snowpack_Layer in July only (“Early July”, “Late July”), stats included, facet by DNA/cDNA, “GL ICE” removed, “Gl ice” samples removed, “30 July_SWS1SE_Slush_R”, “30 July_SWS1SE_Slush_D” removed

plot_richness(subset_samples(ps.final, Month == c("Early July", "Late July") &
                                       Snowpack_Layer != "GL ICE" & 
                                       !(Sample_Name %in% gl_ice_spls) &
                                       Sample_Name != c("30 July_SWS1SE_Slush_R",
                                                        "30 July_SWS1SE_Slush_D")),
              x="Snowpack_Layer", color="Snowpack_Layer",
              measures=c("Observed","Shannon","Simpson")) +
  geom_violin(alpha = .5) +
  geom_jitter(position=position_jitter(0.2)) + 
  facet_wrap(DNA_or_cDNA~variable, scales = "free") +
  # stat_compare_means(comparisons = my_comparisons)+
  stat_compare_means(comparisons = slayer_comparisons,
                     method = "t.test",
                     var.equal="welch",
                     label.x.npc = "right", 
                     label.y.npc = "top",
                     symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1), 
                       symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
  labs(color="Snowpack Layer") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 20,
                                   vjust = 1,
                                   hjust = 1,
                                   size = 12, face = "bold"),
        axis.text.y = element_text(face = "bold"),
        axis.title = element_blank(),
        strip.text = element_text(size = 12, margin = margin()))
## Warning in `==.default`(Month, c("Early July", "Late July")): longer object
## length is not a multiple of shorter object length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length
## Warning in `!=.default`(Sample_Name, c("30 July_SWS1SE_Slush_R", "30
## July_SWS1SE_Slush_D")): longer object length is not a multiple of shorter object
## length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length
## Warning: Computation failed in `stat_signif()`:
## not enough 'x' observations

## Warning: Computation failed in `stat_signif()`:
## not enough 'x' observations

## Warning: Computation failed in `stat_signif()`:
## not enough 'x' observations
## Warning: Computation failed in `stat_signif()`:
## not enough 'y' observations

## Warning: Computation failed in `stat_signif()`:
## not enough 'y' observations

## Warning: Computation failed in `stat_signif()`:
## not enough 'y' observations

# ggsave("july_alpha.svg", width = 12, height = 8)
# ggsave("july_alpha.tiff", width = 12, height = 8)

All measures for Snowpack_Layer in “Early July”, stats included, facet by DNA/cDNA, “GL ICE” removed, “Gl ice” samples removed

plot_richness(subset_samples(ps.final, Snowpack_Layer != "GL ICE" & 
                               Month == "Early July" & !(Sample_Name %in% gl_ice_spls)),
              x="Snowpack_Layer", color="Snowpack_Layer",
              measures=c("Observed","Shannon","Simpson")) +
  geom_violin(alpha = .5) +
  geom_jitter(position=position_jitter(0.2)) + 
  facet_wrap(DNA_or_cDNA~variable, scales = "free") +
  # stat_compare_means(comparisons = my_comparisons)+
  stat_compare_means(comparisons = slayer_comparisons,
                     method = "t.test",
                     var.equal="welch",
                     label.x.npc = "right", 
                     label.y.npc = "top",
                     symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1), 
                       symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
  labs(color="Snowpack Layer") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 20,
                                   vjust = 1,
                                   hjust = 1,
                                   size = 12, face = "bold"),
        axis.text.y = element_text(face = "bold"),
        axis.title = element_blank(),
        strip.text = element_text(size = 12, margin = margin()))
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed

## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed

## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed
## Warning: Computation failed in `stat_signif()`:
## not enough 'y' observations

## Warning: Computation failed in `stat_signif()`:
## not enough 'y' observations

## Warning: Computation failed in `stat_signif()`:
## not enough 'y' observations

# ggsave("july_alpha.svg", width = 12, height = 8)
# ggsave("july_alpha.tiff", width = 12, height = 8)

All measures for Snowpack_Layer in Late July, stats included, facet by DNA/cDNA, “GL ICE” removed, samples “30 July_SWS1SE_Slush_R”, “30 July_SWS1SE_Slush_R” removed, “Gl ice” samples removed

plot_richness(subset_samples(ps.final, Snowpack_Layer != "GL ICE" & 
                               Sample_Name != "30 July_SWS1SE_Slush_R" &
                               Sample_Name != "30 July_SWS1SE_Slush_D" &
                               Month == "Late July" & !(Sample_Name %in% gl_ice_spls)),
              x="Snowpack_Layer", color="Snowpack_Layer",
              measures=c("Observed","Shannon","Simpson")) +
  geom_violin(alpha = .5) +
  geom_jitter(position=position_jitter(0.2)) + 
  facet_wrap(DNA_or_cDNA~variable, scales = "free") +
  # stat_compare_means(comparisons = my_comparisons)+
  stat_compare_means(comparisons = slayer_comparisons,
                     method = "t.test",
                     var.equal="welch",
                     label.x.npc = "right", 
                     label.y.npc = "top",
                     symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1), 
                       symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
  labs(color="Snowpack Layer") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 20,
                                   vjust = 1,
                                   hjust = 1,
                                   size = 12, face = "bold"),
        axis.text.y = element_text(face = "bold"),
        axis.title = element_blank(),
        strip.text = element_text(size = 12, margin = margin()))
## Warning in max(data$density): no non-missing arguments to max; returning -Inf
## Warning: Computation failed in `stat_ydensity()`:
## replacement has 1 row, data has 0
## Warning in max(data$density): no non-missing arguments to max; returning -Inf
## Warning: Computation failed in `stat_ydensity()`:
## replacement has 1 row, data has 0
## Warning in max(data$density): no non-missing arguments to max; returning -Inf
## Warning: Computation failed in `stat_ydensity()`:
## replacement has 1 row, data has 0
## Warning in max(data$density): no non-missing arguments to max; returning -Inf
## Warning: Computation failed in `stat_ydensity()`:
## replacement has 1 row, data has 0
## Warning in max(data$density): no non-missing arguments to max; returning -Inf
## Warning: Computation failed in `stat_ydensity()`:
## replacement has 1 row, data has 0
## Warning in max(data$density): no non-missing arguments to max; returning -Inf
## Warning: Computation failed in `stat_ydensity()`:
## replacement has 1 row, data has 0
## Warning: Computation failed in `stat_signif()`:
## not enough 'y' observations

## Warning: Computation failed in `stat_signif()`:
## not enough 'y' observations

## Warning: Computation failed in `stat_signif()`:
## not enough 'y' observations
## Warning: Computation failed in `stat_signif()`:
## not enough 'x' observations

## Warning: Computation failed in `stat_signif()`:
## not enough 'x' observations

## Warning: Computation failed in `stat_signif()`:
## not enough 'x' observations

# ggsave("july_alpha.svg", width = 12, height = 8)
# ggsave("july_alpha.tiff", width = 12, height = 8)

8 - Alpha diversity with rarefaction to the smallest number of reads in any sample

ps.final.rare = rarefy_even_depth(ps.final, 
                                  sample.size = min(sample_sums(ps.final)),
                                  rngseed = TRUE, replace = TRUE, 
                                  trimOTUs = TRUE, verbose = TRUE)
## `set.seed(TRUE)` was used to initialize repeatable random subsampling.
## Please record this for your records so others can reproduce.
## Try `set.seed(TRUE); .Random.seed` for the full vector
## ...
## 25951OTUs were removed because they are no longer 
## present in any sample after random subsampling
## ...
plot_richness(ps.final.rare,
              x="Sample_ID", color="Sample_ID",
              measures=c("Observed","Shannon","Simpson")) +
  geom_jitter(shape=16, position=position_jitter(0.2)) + 
  facet_wrap(DNA_or_cDNA~variable, scales = "free") +
  labs(color="Sample_ID") +
  guides(color=guide_legend(ncol=2)) +
  theme(axis.text.x = element_text(angle = 20,
                                   vjust = 1,
                                   hjust = 1,
                                   size = 12),
        axis.title = element_blank(),
        strip.text = element_text(size = 12, margin = margin()))

plot_richness(ps.final.rare,
              x="Month", color="Month",
              measures=c("Observed","Shannon","Simpson")) +
  geom_violin() +
  geom_jitter(shape=16, position=position_jitter(0.2)) + 
  # stat_compare_means(comparisons = my_comparisons)+
  stat_compare_means(comparisons = month_comparisons,
                     method = "t.test",
                     var.equal="welch",
                     label.x.npc = "right", 
                     label.y.npc = "top",
                     symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1), 
                       symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
  facet_wrap(DNA_or_cDNA~variable, scales = "free") +
  labs(color="Month") +
  theme(axis.text.x = element_text(angle = 20,
                                   vjust = 1,
                                   hjust = 1,
                                   size = 12),
        axis.title = element_blank(),
        strip.text = element_text(size = 12, margin = margin()))

plot_richness(ps.final.rare,
              x="Fox_Aspect", color="Fox_Aspect",
              measures=c("Observed","Shannon","Simpson")) +
  geom_violin() +
  geom_jitter(shape=16, position=position_jitter(0.2)) + 
  # stat_compare_means(comparisons = my_comparisons)+
  stat_compare_means(comparisons = site_comparisons,
                     method = "t.test",
                     var.equal="welch",
                     label.x.npc = "right", 
                     label.y.npc = "top",
                     symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1), 
                       symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
  facet_wrap(DNA_or_cDNA~variable, scales = "free") +
  labs(color="Fox Aspect") +
  theme(axis.text.x = element_text(angle = 20,
                                   vjust = 1,
                                   hjust = 1,
                                   size = 12),
        axis.title = element_blank(),
        strip.text = element_text(size = 12, margin = margin()))

#now without stats, only shapes
plot_richness(ps.final.rare,
              x="Fox_Aspect", color="Fox_Aspect",
              measures=c("Observed","Shannon","Simpson")) +
  geom_violin() +
  geom_jitter(aes(shape=Month), size=5, alpha=0.7,
              position=position_jitter(0.2)) + 
  facet_wrap(DNA_or_cDNA~variable, scales = "free") +
  labs(color="Fox Aspect") +
  scale_shape_manual(values=nmds_shapes) +
  theme(axis.text.x = element_text(angle = 20,
                                   vjust = 1,
                                   hjust = 1,
                                   size = 12),
        axis.title = element_blank(),
        strip.text = element_text(size = 12, margin = margin()))

plot_richness(ps.final.rare,
              x="Snowpack_Layer", color="Snowpack_Layer",
              measures=c("Observed","Shannon","Simpson")) +
  geom_violin() +
  geom_jitter(position=position_jitter(0.2)) + 
  # stat_compare_means(comparisons = my_comparisons)+
  stat_compare_means(comparisons = slayer_comparisons,
                     method = "t.test",
                     var.equal="welch",
                     label.x.npc = "right", 
                     label.y.npc = "top",
                     symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1), 
                       symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
  facet_wrap(DNA_or_cDNA~variable, scales = "free") +
  labs(color="Snowpack Layer") +
  theme(axis.text.x = element_text(angle = 20,
                                   vjust = 1,
                                   hjust = 1,
                                   size = 12),
        axis.title = element_blank(),
        strip.text = element_text(size = 12, margin = margin()))

#now without stats, only shapes
plot_richness(ps.final.rare,
              x="Snowpack_Layer", color="Snowpack_Layer",
              measures=c("Observed","Shannon","Simpson")) +
  geom_violin() +
  geom_jitter(aes(shape=Month), size=5, alpha=0.7,
              position=position_jitter(0.2)) + 
  facet_wrap(DNA_or_cDNA~variable, scales = "free") +
  labs(color="Snowpack Layer") +
  scale_shape_manual(values=nmds_shapes) +
  theme(axis.text.x = element_text(angle = 20,
                                   vjust = 1,
                                   hjust = 1,
                                   size = 12),
        axis.title = element_blank(),
        strip.text = element_text(size = 12, margin = margin()))

plot_richness(ps.final.rare,
              x="Snowpack_Layer", color="Snowpack_Layer",
              measures=c("Observed","Shannon","Simpson")) +
  geom_violin() +
  geom_jitter(aes(shape=Fox_Aspect), size=5, alpha=0.7,
              position=position_jitter(0.2)) + 
  facet_wrap(DNA_or_cDNA~variable, scales = "free") +
  labs(color="Snowpack Layer") +
  scale_shape_manual(values=nmds_shapes) +
  theme(axis.text.x = element_text(angle = 20,
                                   vjust = 1,
                                   hjust = 1,
                                   size = 12),
        axis.title = element_blank(),
        strip.text = element_text(size = 12, margin = margin()))